{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors\n",
    "import matplotlib.path as mpath\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import datetime as dt\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_data(fili, variable, swap_latlon=True):\n",
    "    \n",
    "    ds = xr.open_dataset(fili)\n",
    "\n",
    "    if swap_latlon: \n",
    "        ds.rename({'lat': 'lonx', 'lon': 'latx'}, inplace=True)\n",
    "        ds.rename({'latx': 'lat', 'lonx': 'lon'}, inplace=True)\n",
    "        \n",
    "    da = ds[variable]\n",
    "    da = da * 1e3 # Scale to mm\n",
    "    \n",
    "    # Get attributes\n",
    "    da.attrs = ds[variable].attrs\n",
    "    da.attrs['units'] = 'mm'\n",
    "\n",
    "    proj = ccrs.PlateCarree()\n",
    "    extent = [da.coords['lon'].values[0],\n",
    "              da.coords['lon'].values[-1],\n",
    "              da.coords['lat'].values[0],\n",
    "              da.coords['lat'].values[-1]]\n",
    "    origin = 'upper'\n",
    "    \n",
    "    return da, proj, extent, origin\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_latlon():\n",
    "    \n",
    "    diri = r'C:\\Users\\apbarret\\Documents\\data\\SnowOnSeaIce'\n",
    "    latfile = 'Na25.lat.361x361x1.float'\n",
    "    lonfile = 'Na25.lon.361x361x1.float'\n",
    "    \n",
    "    lat = np.fromfile(os.path.join(diri,latfile), dtype=float32).reshape(361,361)\n",
    "    lon = np.fromfile(os.path.join(diri,lonfile), dtype=float32).reshape(361,361)\n",
    "    \n",
    "    return lat, lon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "diri = r'C:\\Users\\apbarret\\Documents\\data\\SnowOnSeaIce'\n",
    "fili = 'era_interim.PRECTOT.201501.day.nc'\n",
    "\n",
    "precip, pproj, pextent, porigin = get_data(os.path.join(diri,fili), 'PRECTOT')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x10f53da0>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD8CAYAAAB+fLH0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztfVuofttV32+sb//PJRdNUvUQNdRY4oMtGiWklkqxFWua\nh0Zfgj7YFAPxobUKffBYaWsRwRYvFIpCRPG0eGlAJUFsQwwWEawa0xhjYkyMCRpOcry0mERPzj57\njT7MNdcac64x5hxzrvXtvf5/vwGbb33zvr695m/9xmXOScyMi1zkIheJMtz1AC5ykYscSy6gcJGL\nXCSRCyhc5CIXSeQCChe5yEUSuYDCRS5ykUQuoHCRi1wkkbOBAhG9iojeT0QfJKLHz9XPRS5ykX2F\nzhGnQEQnAL8P4GsA/DGA3wTwjcz83t07u8hFLrKrnIspvBLAB5n5Q8z8DICfAfCaM/V1kYtcZEe5\nOlO7nwfgj8T3Pwbwd63CD109hx+995ntvagkp5H5VIsXCvTWNeu1lheZtbEIRmgXdbYF4Iu+5C9F\nrXUFqwlv2STNMXY7naplvPW9ddS6SuF1ElXyK31Uyn3oPX/5p8z82bWy5wKFqhDRGwC8AQAeufcZ\n+Ht/65vTApZaM3K9XJ6Wfae8jbxMpX4yhlLZcVTTVypb0p5eZy4j87O6rPWtpLHWrmibHeN561vf\nNV/fTGVG8RiPGEU+r9OmsqNo8yapv8hNUkbLF7cj08WEuWESZdKJJMuNbOel5QazXPie5svyedkw\nvqy8QuK7601lXvuyd35klanIudSHjwJ4ifj++VPaLMz8RmZ+BTO/4qHTcwDKbphonaaJVkZrS/Y9\nEHgotF2pD1m3NMZhUMtRsb3Kv6SWvyrve5OkXbTXuQ25qRfpaFMHhBuQCQilNsJ3GxDU+pWJnY8l\n1mkBhBY5Fyj8JoCXEdFLieghAN8A4C27tKw9sL3g0dqGR0rtlICh1taOE9XV9xnlxkWMp7IHXLBX\nmuQ5IKzz15P7XNLb9llGxMzPAviXAN4K4H0A3sTMv1utqLED7QH2AIODQVTZgizfwhaSyXyGn1iy\nhRawOMdYziBjvUi1nkd1KKkNUlrUhmQ8PJQBRHnTayyhq94GsDmbTYGZfxHAL3ZVJkp14TjJZNpA\na/uCp15eJq+fly+1780zyhHRot/L+6FhbTtwSNLeQcT7pm9hD5Y9YS/xqg3reiUG0c7OPKzCqzK0\n9H9/vDos2YFSF9kC0K9W7MUWzk31b8l+MHa//48jrSyhJB47Qo/soY4cAxQsqt+jSnjq1dQIp10A\nA/mNjntIi0GyAkR3bVdoFc3rYJftVx28k7FmXCyW7QCELSpDK8AcAxQAsPWQerwS+eS06hW+b/JG\neOp5PBG9b+1zsZ0zy1hRLaLqYXkdNFfkXoCQW/y9xsXcjtAKCGr7nV4GzWvhkcOAAhCAwQSHPaSl\n7T2AYIu0uh7/Gsu5FZMWtWGr9Hgn9vZgHPLJWwHDFlWiJLehRmxhCxEYavnJsDY8sBoQVcDppsMg\nWpL5jV9lETv0VfA4mP0W1IZqQFOj2tCrMmgMoWbjkHJIUAAMdcKrStTqFNooAsNd0vBa317VQ7M1\n7GxsjEZFbWJrXgav56EFfjyqQzoGS70oqAKddgSP67GXIWj31wIIwIFBASjYGWrSCgx5v95J4mUS\nrWyhJpXxzW3dJ3EJUbZyjnPS+nVf/qjF2gTvcT1udTuW5DhPjTEStyqxB2OwJJ/g5oTf6aE8Z9jz\nfSa10OYakPQYGJcy5fUKWrlq2Q7XYy8g1IKnLDnWEzVAHZFLldDSNgDDan1Ei31BS3ewhar0MKfY\n10E9ED1SsyfcKLaCPQHBa0fI1yd4VIaq3eFMKoOUY4FCFC8wuNo602QoLoRyAIPa5MaxNt7rXcYq\njLO7UbExKLaIGiNQJ8YZ1Im9oxZLsicg1NZkSDkmKAA+YPB6JUr5lTUOXWzBIzXbQk2FOBvY7f9I\n7BnN2BKboI7ljGpDOk77jW+terTKa/nW2LYCAnAkUFBdjDu1tdFVublcr7gNnju7Jh8w2esN3uJ+\nlFILUDoSIABHAgXABQxqgJOHMWzwSKzYgsUszmlbqMUseKTBK3HbeypoXKIU2txiYKxJD0uw6vfY\nEaTcNSAARwMFwAYGBRyqdbcAQy1+oVRP66/hzV19y6v3Pqz7vAMpqQrabkt2WSXtDAZGDyDcYDBZ\nwhZPwxEBATgiKADuCeQyPtaAobPvTaqDxRY02WOSe+IVDrLbUuumKt63ec+iJ72dvniEWp/nBoSa\nrUXKMUEBMGIP1sU2qxJaeaPurmxBAwavCiFFG5Nngh8gsMkTybhFdWiZ/CWWYImXJdRcjz2AoMUg\n5Ixmaa8NBO/+yYhiTU4NGLx2hmI7FdXByNvVvqDJlL+HofB+MzZu2XotPviaLcHjcfCICRyVSZ6O\ns01l8K5jsMCgFRCAI4FClHM+yC2qhBcY9hzPKjsaBh1gCXS5E/cCjjF767dMcNceCbwuO9sFOh58\nsx+DJXiDgbZsltK72nEPdiDleKAAOMOYUTc+tsYxVGIWkr48akGLp8BrW5AT3zGhV5P+jthD02at\nrjIKI6iwAI0lSBruAYTVLs9OT0MtWtHLEDSVIRcNEDz7NkQ5BCiooOad0FtVCW2DFqusN3ipBgy1\nNpsDoRx2hQ4bwjnckrnnocQSlmXUW/pb2xd2i1twGhZrS6A1QMjFYz+w1IUWQAAOAgpAQxjzXkFO\nO0QnblIj1PtwBiHdR4ug9oxmVGMZCm7Ilv0SvCxh6bdsG3D361AZtqxj6BnXpqeLiD5MRL9DRO8i\nondMaS8iorcR0Qemzxd623O95a20rapEzfAYv5fiF3rVCK8nogRCmnrhlLswSObqRG6DaIUSr8fB\nu8WaJ3/p26derMruZFTUGIK2Z4NX9njl/ENmfjkzv2L6/jiAtzPzywC8ffreJP3LpdftbFIlelyV\nXnXDU9bjiWiZ0DkTsVSKO2QiuT1BUx1WE6fB47DkrxmApOSWHaHkaWjZWKWmMnhdjl51ocWIeY7/\n/msAPDFdPwHg63oa6V8u7Wyr1m5juerCqb3fxqWJm+XtZmw8gIuzZfFTqW4pbQ8p2hs61IEtBsXW\n/raCAgP4JSL6renAWAB4jJmfnK4/BuCx3pGYwFCNQVi3x/lb2atKlMp5JolWXjKSVttCafXkyri4\ng8HyDGJ5Ikpeh1ydqIU052k1taHVjiDf/OdWGXoBQfN4eGTrCVFfycwfJaLPAfA2Ivo9mcnMTETq\nE5CcOv3QdAx9HL886DhS6NzvTeQ4RWrdFnlPi0pObMralvXENQ+0nGhtlFH7jGmlPGtstfvYQzpP\nrJr3ajRAwGoxTy97KGzVwas2WOJ1PSZ1GiIVe+0Huex1sGyUTUyBmT86fT4F4OcBvBLAx4noxQAw\nfT5l1J1Pnb539ZzqqPo3WWlopyW46VwyM4KdNDvrHg6gCkSxAp1ye4Jr6zIng0jzy0FKva5Hbxtm\nv3cACMAGUCCi5xLR8+M1gH8M4D0Ip0u/bir2OgBv3mtkexogzfIlYCh5F6bvyTZuudcil5odIvNE\nFN2Tq9/B+Nfe4ZoHz+pIj3hUhxbjYil/1UbBsGiWq6gMPQZFy7ugr5vQ10RYskV9eAzAz08P6hWA\nn2Lm/0lEvwngTUT0egAfAfBaV2sqTZ4+a+qEpjrk7bWoEqu6mSrBDrWidI+ynueg3NqYpBjptcNn\n53yr3TNJjSWUbAm5G7LEEmoT3mtHaCmT52l91+wHW9hB7/LpblBg5g8B+FIl/c8AfHVXo9rktvrP\nJ7XVXgswFOsawGCNz7IvaBLbdrS7+8nSwwCM2wKM4rqHLYFKLWdF7iUtk6bH578VEHr72dIWcKCI\nxkQc3gSgd9l0oY29VInYdi1QqeS6LNgWzIVSVvpQUD/OIOrGq9lGra2T34pN8HgcvCzBikfwBie1\nxiGk+X0BSV5AaInuPAgoGJPD6WpsXzadtlOs7zVM1oChxb5g9W2NxbIf3Efh0IC91sGzGjL3OJTB\nQg9Smss5VYY9AKHXfqC5Gy3bQevxeMd5aswHXptA6yQXMJSCj7yTf+9j5Tx1PMbBFgDIQ6kPsOGK\nJrbLMp9k2xnQpt2YNjCE1TgUQFiX8bODnt/mWE+DNXEdHgWgZ73Duj7Lt7k14T2uSkuN0NhCTvtL\nnoi5ilOF6GULZ2QZ+SGynqXSS13xtq54HLRl0S3ux1bDYgkQasueLXUhH9u52IGUY4FClBI4SNlD\nnTDaUMtawJCrBsrEr9oX8nYtVUN87zl/stWu0Lt8uha4VK+Pqb6uOlgeB/VMhQY7wtL/voCQ5g0r\nQEjGttF2oHth/FP9GKBgPXceYAD61YlKG9Ux9TIG7315xpCkGzELmhxEZci9Dr17J3jof8mOIGVP\nQLDGsNR1jHujMfG+PQyGqbDZipbmZA3FtgrA0KVKeBhGntbCFqyJbIHTHRwjV9qGrZc15PsvtjKB\n8KlNLH3C1/J7VYZeg2La/lpdKLGD+/8wGBjgsMHWsFIntAlesTMkZbVra4GTMuaifaEFPChTB4q7\nRw31MlJUVaf/UfFGMuYBS9Vt0U2VIgWEtI5tR8gXOeX5sYyWntc5l/1gfT8+MGjZqOVwoBDFzRo0\n8dxVqa0eVcJKL034Ul3z7V8ZXK5C7HXORYeUmEHvzs19uyrZLKC2TmBvQMjlHOxgqxwHFJSH180a\nHOpE1QCpMYZBqdujSmRpq/URVnnTs9BgO1BEVRM0sNnRA1EKWtKPjIuTNlUdzHIVlqB6Ipybpazy\nnYCgGRTTMO3zswPN61GT44ACsN75aJLN4JC0ValTsTOo5WrA0GpfqNk+VtUzVqBNZvUgmXXabUQ9\naq5Ia61DSMsnfvpdi1wM6UpgUGHCuyIaGwAhKZeBwV7swHZJ9k/tY4FCFIPyNhkik/bydhTWUGuz\nVq4FGOI4PG7KksdCpLsns9c1ektyjrUOrn0KlAmviQUItZWOaV/l33cLO1iX05lBS3DWMUEB2Ic1\nJO3Bb0TMv1uqRF6u5q5UVBDTTVmzVwyDw75gM4gjnh5Vi00A6izB422w1jTM15WgJm/YsqYu5CBz\nU1A34v1uZQet0ZqHAQXzzb2FNTjUiWL5wsQpAkptXDWpqRFFI2lBhbiPpOaG9ErJE6H3u7YztACC\nOY6OqeY1ppbYQc/vtnU7tl3F3HotPujZOv/4myUbvsUJ49xbwdyfQdsjQdSLdSkvJ/uX45b5WZ15\nmXU+zsr95+NV90Sw2jyjyCXU0h2ZHvnmF68tIWcJ6slPGUuouR4tQGjZB8FjP5DS4lXYqipocihQ\niCLfwpTvaQD4waG0EUv8LS1wkOWTib6utwKG/DrfLyFvWwMGDURku/Ez7odQm/wSKJLrqf4ZRbom\nbxSQuGF/bILHuKilh372B4T7AQzuX++D8SyoaoVXpfAYIT39aeLxSmjiad9SEzxuzKSdxn/vhtDn\nfIOVYoxCpa2aG1JvUy/rjUeI0soQkr4KdH/rvolSzgkIwNGYQryv7HlSKb6XNVjqREzL3vyxP/Pt\nL9vK1JCqKmFtuybYAoBlxyZtvJqIt31RhRgIGAeAR2Ag0IhlF6eViqWMt1Es1UGKXOuQTD7HVmt5\nPUttUFlAKa+RIXiiE9PxlhnCXbADKccChSg94AAkDzFTo62hpE7kdQ07Q1WVkGpEbM+rJsj71dSM\nCAzStgDAtT17pkL0bvnWG6VYE2lL8AQpqemKp+HcgLAHGLR4FIos5r5zSVrjbVEp1HJam1T+XvNO\naNctqoQjTqDqpjzn1vMbV09692kcjevWU5/mesY/znWi8wYbgtaO3m/dfuCVVkNia/njMAU5AZI3\npigjkr2sweWhqKgFST+WOiGYhluVqBlDZVqeno9FSzN3oXawB9nkQOAddnjOmYRHdSgd7LLULasN\n51YZWtQFjzFxD3awxQNRfS0Q0Y8T0VNE9B6RZp4sTUTfSUQfJKL3E9HXdo1Kiy8AAkBkyTGYaMUe\nsrepGfRk9Tugcyn1Ul7Nz9/2htEwWR+hjTnf+5EoWRORhD7LRVJZm8UgJlm3U1anS8PHEhLPQAIS\nwScfVYlzAUK6VqIckLSM3w5G8gQh6bszEaTaVCory+ey9yrJnwDwqixNPVmaiL4YwDcA+NtTnR8m\nopN7NLlYD6xXrVCiIt3Lsuc2ykNM2snquFWJHBgstSH/3KJG0HBeNUTIbW7bnvZbp/qlwKQt7ddW\nK3rYQasR8dZcksz8KwD+PEu2TpZ+DYCfYeZPM/MfAvggwlFyVSlusqJNXIU1hHZ2YA15f4I19DAG\n1ia1c+IXd2vSRNgE2rZrG/x9dEh585Xl7TsmxkRKWIIsHz7XLEFfaWiwh4rKYDGEUrjy+u2fMoQW\ndpBLD5Po8UL0ckPrZOnPA/BHotwfT2luMcEBKIODoVbMElmDAIhVX9ok9xoUtQmvlc/raRu0lIBB\nYws5cIoJTrEPqQa0qBAFueERI1iNURiZy7EKU8CSZksAkE2khdJ7ohZDW0MyWSxAyD0aEhCWsfrV\nBTn+EhjIscv79E7uFuCQ9+uRzd4HDr6rZisUEb2BiN5BRO945vpT63Zr4KCmK+NzPPCuoKcoFjC0\nlC/V87TpaUetd142UJPcnrBfu/rEn/PzCVWZIF5AkOVL6oK1z4F2D3k5vb1GteKWIhqtk6U/CuAl\notznT2krkadOP3TvuaZuH8HBtAU41IrNrKFHnagxjJp9ILu/ohrh2SdSq6PZFSaD5TlXUa4Oe8lY\ngRWXYLEEKZZhMVJ+j1ExjklTF0rsILQlx7+NHWjlWtQETZXySC8oWCdLvwXANxDRw0T0UgAvA/Ab\nTS0XDH9Vu0OShhVAmPaGGjhY/ZTUibx8BiSrcnIshr1hpUYk7WeAlHsiNhoV+7d6D/aEPLRZUx00\nuwGguwU1O0IJEOR3PQgqtR8s4/KxA01dkONfb+K63o9hLzWhFwyiVOMUiOinAXwVgM8ioj8G8O8B\nfB+Uk6WZ+XeJ6E0A3gvgWQD/gplbFsTJjsOnYqBaRSvKOppBi7CKcVBXYpZWIcqxyH4G+CMak3ug\nJVrSE9Fo3Vur0DTgLL6Dxg4dUJHWHZuL6wUMkPCkeVSGpq3eK+zAytPE4+UolelRE1rWkVRBgZm/\n0chST5Zm5u8F8L3uEUQRkysR+dYUD7J8ltTApKx8HjpdW0+htm9N9Pi/GAvh0bG8BiL5xLdCmYkQ\n/+80inZjHwOFtmN6vgKytp6htmKSBqCA8bVoxsAa1v/m3MBYUhtq8Qhpu7aXwWtQbAlEKp0+5bUZ\nNIcv7wQEUo4T0QikysyeAOEAh1BXiQQU7RNjPdFl+zXWMFN9nu+xGRjy36O6bDtMdALS9QxE03gH\nM7qRiMDON/8qWhFchIjohtSiF6OY0YSFAKWSypC3qwGChx20gEHet7eMVS4f16rOhkjGKJu9D2cT\noYOrUrA7uMrWfjtjOzi1zZrtwBqDx8Zg9KueH2HVaZF87UNHe9o5D16vQz45TNuCEbGY1C2oDJZB\nUau7hR3UAKE1ArEXEFpsDMdgCpMBcKXnA2X2YNgdiusdZHn5G3KBNQC6SuGwNazUiRXlh4MxyHEu\neTzZAkw1IvttKKofOTugARjGxS2Qqx6NayWk5EbG3MAILCwhD2de8tcsIbS1ViUslUFrr2ZQ9NoO\netjB0YBAyjFAYZLcM7A2Bopr+YyeUa3QwCG2re7ZkE/4aazJpJdlM2CYh2KpElJKwCDvVSyrXtKn\nToflnlxLpjcYPJMFUJPqoAOE9BKkdoQelaEECK3soORi3FNN6AGCvTZyOa76gEpwkDVyw6XpWhAF\nrNSKUsh0NVR6NQYlX6odWlsOVUK9Nlybe7gnb0N6V/mVJsaegKC5GKX0xiOU7mEEdQPCrt6HW5Ps\nTR2lyB7y36GBPaze8rJcA2soqhQlw2Ls22AMRVVC8UiQtVHrPIaCZ2EgRI4fjIv9Eu0JMT5huQ75\ncp3DUkcPzIksQVUPGlSGkkGxBwzSMe6jJpSAwJJzGRyPAwpRMj0/F3NTV6CuXpTAQSu3VaXIJ380\nLI6ZnWEelFImppeAAdDVCKEahO/CExFtCzFmIdoVRg7lbm6Wz4osJz6ln+F6kXTJNK2iF3O1Ie1j\n/VZO0gyVoZUd9IBBt8fhQEAg5dDqgxayLKVZvVBUC98+C+k43CpFizphlEnSLfWj1I9cMHVQKcXy\n5zYDzR2puR3D97VBcc5TIhk9gGB5C0r5ewBCKUqxpFbU6mpyCKbAUN7CUgrswa1eOJlDaAPrN7gc\nh+apUFQKrzoxN62U4VFRJSIDyNjIrEYYqpNMnz0RNzFdj4GYYxUoFvZLLC1XREoDo/Q4LHVS42IU\n7Q2vAUKNIdTUhV52sEVN6LET1FjBWcOcb1uqHoj8tyiARBUgWu0OHrViAzjMo8nsDKoqITwbyf1L\n+0Ic+6Ql5Bu8hvwBcXdn3GD6bLcqxGjGWYUQ9oRcLI+DBATLjlADhJK7UXM11sCg1WZwV0Cw5xby\nx1YfAH0Rk5SKeqGfG6G1o9H7Spns62olptZOYTxmEJM2hsJvwpbKUFFn4p4LqxWSBW+Fd6PWeeKK\nH61GxfP6HoawHl8dEJLypUnZamRUaLtF9XvVA49qoG0FV5LjMIWCigBUGEQPe3Ayh6pakY27iTXI\ndjTPQ5SqKoGV8ADQjWJ0FHELSTBToBKiz4K3oiKj8Zm8vZPrNUvQ7Ajz94LKYKkLe7EDDxBoYgGB\nt6yn3jKmbfajYzIFQpEBABUjY6G+yRzyX8IwSq7K5P1q/ZRiG2Q7OWNQ8tWNWqw3ecMbf6nT/0ho\n4c2aJqJ5HEL62o6QqxE1G0Lat0LbK+xgb0DQ3vI1VmBJjRW0MgJLjsMUNGlgAEDF/iCyttgd3MzB\nwxpi4TxWASWXZJYX8/NwaET7gkyQzEHYFqQ9YeTUriDckjRQaZFkuL3pUxbLDYwxLf1M7Qa5YTEH\nBK9BseRZsCZ8s/tRURFqZUpla3WWsXhVrzagPwwo5Pen75eQV8rbcKoYWwEiUy30cGfZnehjNhAq\n4DBXyNQJb16uSsgwaK3MMICYQ2hzXBYtApm8cQojRmiHyMpl0uP8Zl8CleTElmpDjSGUAMGjLnjA\noJcV1Mpo5WrllzHUQaAVADQ5pvoAI34gl4qaYRoajXpuw2Qt1sETE6CpFLlqoBkg45/Mq8U3VIQU\nlaTF2KiJx+sQ0lJGkOQ5AUHbag1I6fRegOA1Hu4JCB61oLZPZIscgykUJ3ZWdCOD8BgoXYbJglFy\nbUwU7s4Ca1CNkMLIOA95NjLymjHMbYreczVCsoVZLREMQcYkEPkWSgHId12SsQlSXZAuyHzTlHkX\nZgMQWtWFXjDQlmmn9+ZTEXpBoCYeAOi1LxwDFAD9TaTsEqTd5wooCiDRqmL0AER5PcTS/mxzmO0I\nvK4T65kqg/BMJOmoA8PNErdAADiunEQAiZ4t2uJSaWlgzFUH6W0AUjuCFxBawaAVCGTb+X20ltHK\nJXUcLMAjexgZgQOrDwDqG51MUlU1CmpGi4qhb/ya10vpvB32nLYb2iK9jqYy5H1Fz4TmuZjSVsfQ\nSaNjck/7PFxy8ZN2sItlR7BUBiB1ay791NlBlB5A8HoQWgGhphZ4VILYxl6AAByIKZibsQLdLGIX\nBuFhDxWjZMIcDItnVCuqrEFTJ0RewhikKqG6MjEhalRHKMQsSBfDMIT8zOAYD4IJzCAc/hLUhdBd\nckak8DqkuyevGUJoe8iAo8wQ9mIHORDkskWNWPorg0BJWib+FpA4DCgATtUgigMoqu05QUIFiKl8\nFSAM1cJSK3rAoQgMwDL5gbU3QrAFYg7BTEzId2Miokm98K9/iG7InCXM3gZDZcgNiiXvQmjbthtU\nAaLytt/kbtygFtwWAGhyKFDQxNw9SRPlKHqrvaLB0jBUqnEQCpCsbA+5UQ8ZMxLpsc05snGKHUgM\nkfmyaaKVO3UGBmXvhtm+kC/F9kh1R+eAIzcK5Y6TUa5tkLLyPhQAoYUd9DCDrazAmqh7GQhbgaDF\nM1EtaRxF/91E9FEietf092qR13cUvdSdDYm2g/xPvzNK/xxtLWPJ/uY6lPxZZVdlBvEn7jPpW97/\n1NbcxnQPatkTzfYEHgbM50hKG0NsQ/Y96HlEkSWET6r8TwAgnhupHSKbexuu+SphBs/wCTc8LOnC\nqBjVhWu+mhjEgGs+zewgMo7FcLkcVW+mib9oK5h3lBZ5AJJ8rcxcTuj11hmSln2gVLeljNVfj6vS\nwxR+AsB/AfBfs/QfYubvlwnZUfSfC+CXiOiLmg+EsR5C423mUjs61I3lTZ53GMunGVRYCzHnJ+oF\niaLLWBbyPzGBXKWYxp4oCScE9SSxKXCoOxJoHAOA3Cy6A58GUPxOtEQ5DhP6SG1hGJa64p7kgqjo\ndQhsIV0eHdWGHBA0lSGmA5jBINzyWlXwqAklFaBVfVjabmcCtbd7y9t/r5gETaotG0fRW9J9FL1L\nHGwiSpVJAEUmIdtYjwOqN2O9ZiEttz7TUpYVAEFZmmAO+fjyuvM48t8p83Cs8iJbyM+RHIzrTHLD\norQlxDf/rD5kEYtLvfSttjcg5F6EWp4VDKVt8156I9fe7l7vQe+b/zZXSX4rEf0zAO8A8K+Z+f8i\nHDv/v0UZ8yh6InoDgDcAwMOPvMAXpJQ2oKcrbML6PZI+CkzCZCIKg2hmD4MYhMIccm8FEwWqcAp2\nhIQ1EIWBScYwpfG0mUpgDEO4N6LAFuLqycTVOQAYJ1ViYhCGxKAl6XVY2IGtPsRrqS6EdoYZEFrA\nQLMXbGUF2mQ6NxvoYQFHcEn+CIDvQXg8vwfADwD45pYGmPmNAN4IAM//zM9fzeTaPZqgsREseoHC\nDRKZFTONcFwDBNMy9pJaoYEDM8K8Bhb1AsOiSgDAOAbDI9NidBwmz8O0pNrazDVRHZhnA+M1DwlL\nkEAwA4NQGaT9QIIBEI2K65iFPYBAYwKalBiAJeeY/Ht7GSzpAgVm/ni8JqIfBfAL01f3UfRbpckr\nAaRgUbBo7Y7UAAAgAElEQVS0F70TQLLasVovncNTuWkiF7Z4S2wP0a3JmgdiqRc9FTxyBmxY2Rnm\n2IRoOwCE0TH0Oe+1MK+g1B/g/Hi4m3kyL+wgpAuXY+KF0AFBsoMaM7C8CJ70OV9RCTTZygRuEwhK\nG8aUpAsUiOjFzPzk9PXrAUTPxFsA/BQR/SCCodF3FL2ho89SmfTNrMLBJqoqh8EiVPax6qfMIHT2\nEN7iM3MgWtgJL2rFwhqmtz/zKnaBTxDGx6BKMARQJaHPNC+ljmdGaGbjJWgpVRuu+QrP8AnXOKkq\nwzWfQn0W1wV24AWCkmrQohJYk3LPyd878XsnfU16j6L/KiJ6OcLz92EA3wIAux5FnwzCUaYAHC6b\nAtANFq1AofWSAEEjQEi1IgGHMYKHVCciGGBRJYZwY7N9gRk4nUIUIw2Ta5KDy3N1e8EVOS9+mtSH\naz7NAHCNU6IyhGudHeTMwGIF0kU4j6VBNbDcg5psiUL0tGHX2WfSt2793nsU/Y8VyvcdRT9R4W5R\n9PmauI2bFdWjqMoo6oapamSuzsRIKdUEYFnwRMvJ0IRFrWBQEhU5A0dUH/LrOX4BoGjwjLYFeR+K\nRFfkDSvRi/lkZltd8ABCjRW0MoIjAcEeILDH2Q/HimhsuZ+qd6K9jdL/UD1RatWerj4AAGkBVHnU\nJaVAMXsa5jxe1A4GcIqUXzKHZbwsx3wzzqwhWgEIQ7iWNodJ3cCJg+Fy8lAgbgU/hTon6x4Q8Oka\nA56ZJvc1X81M4enxIVzzabYfhPQ1GFhqQi8QeNnAuaMPQ7ltE36vg148cixQaBHvb1QCjwY7xibA\nMMDCAxSyhASJqGIsLkodIMLYOEzsaSw0PeyMEaATaGBw9FJMOzEFlsLAeAqAMvBibCRJ2UVYM9Os\nNjzN94Itga9wzackUCkChAQDCwg09cALAl4A2GshUs/EP8dk11aStshhQMHaiFU9HKZFdAV+13pV\nNcRQP1S1QzmKLikzE4W1isGIXggkqgUJ3YRBwT05RsPhZIScwqSDjWEAnzhEO0ZDozbBJsYhIxjn\nqEUZk4DBBQiWvaAFDPYAAt/ag743/54gsHXyW3IYULCkuGuzIi4Q6WUZtXqzTaDUtZ7JeWenOOmz\n5ueDZCcVQaoT0i4j7BY8TgnMi1oxAULYBn6KX8AYgAGYYhoA4hPoJoQn0WkIayaE6fgZZjzDA64x\n4Gm+wtPjPYwYEpXh6fHe7FlIvQ8LGMiJn6sGXjaw5+RvnfR7TPZzTfJWOQ4odBgKNWkFEaAAJK3g\n4Shv3lZ0MeblCalx87SUm1NH1kEifuYAEfPEORNEQYWIQQcBGIZQ/eYK9Oyz4LhAaiDEzVqvGRMg\nnPA038PTfG9WH66nxU5P8z0R0XgqAoFcZATYLGCLN2HJd9gMNkz2c07yc7kjgSOBQi7Ot/IeUgKS\nbuZhejLscqYXI1M/VuXijk3CHhGLMGgxUk6qAzECIAwAjYvxkjBMyDFMeUJ9OJ2mdSJykgYjY3RB\nJqsdmfAMX01x91FtSJmBBwz2AoJzAsA5Jv85J31NjgsKNdnr/1ANjOpUXzaOj5dZnXWQqhVLOV5C\nl2W9MSvI098UmxA3UQExmKKBkcLf5JVgZhDdAFdXswcieh4CQwjM4BPjo/jU+NCkPtybGUIEDA0I\n0nBmHQD8bsSCqtDwD9ljkt/lpNbkthZE7SZM4a3VI8kuQl0NbKwPZG/7xgapwkYS94PsJ0tPlmiL\nfGlyiPgwinLDtLCKCXQz2RGGEXQi4GYIag0QbA3P3gD3rsDPPosRI55mxqf4agKDh/GJ8ZHJfnCF\np/kK1+NVAgQRBCQL8BoMrUnmmew9k/y2J/VtrWvwyCFAYYv0/O82A8mqQUeZYsSlU32x1BRaJy11\nKGUVzLMrNAQ1ZfEQAKL7gokwMAPjZHy8OiFuvDLbE/iET40PB3vCeG+2JVyPwRWZg0F8+Os2A0U1\nqPzQLZP/3JP+SJO8Ve57UOiR23wJzADU+YwkgKEBS9buaql2rMZLhRlnhNEyeCYouCpPFKIaGeAT\nYbgaQM+GbU/o08+An30Wn+YRn+ArfOzZF+BPnn0+rvkKn7h5ZAaC6/E0exxydhA+7dWLQH2Cbw4G\nug8n7W0FMB0HFO6//9FaNO/BLQOQS33J3Jwx6CmAx6JvEAM0DhhvGDQyhodOuELgER+7OeEDz3wO\nPvTMZ+OTN4/gejzh0+MSqbhsqBLtBz51YMk/o+X+gXjYzifHAYUHQe7iWZMTvFOVWgGJcHkSA7gK\nrGE4EYBHcQXMgPDUM5+BT49XGGe7gcEEWqj9ZdLeqRwEFIwDWS5SF8/PVigzMfxwLctFwhAXShHA\npxOuH30I189/Lv7tu/8pXvS8v8RnPPw0HjldYyCe/wBgyGjToAVhOGTY3QB0kZocBBQuMsuO2Fh9\nOUtzRQIINO0aPXmGTgEcbh4Crp9PuH4e8PQnH5437ry+d8KJRjx0usEAxkDjAg4GSFhpST7xrXkB\nLuCzyF8fUDgoEdmkOlfqrtpWN4cV3+P1gGnL+HA93gufN48EQLh5lMFPn/A0AjA8+vAzuHcKoHDv\ndIOBGKcJGOY/RJBYJp9kDxqTqIGGt0yx/tTv0eIKtsoWkDsOKBx00kY5i92rsU1zDIrqtQaELD0D\ng6AmLDEjPAAYwiLJ8R6BT8DNI8D4EGM8AXRD4OsBz3z6HgDg+uoG16cBD10NOBHjahhxGsYZEGhW\nLU6LmsHppM4BQ7ND5OAh7Q9dKsqOkbEl6VWfemULyB0HFO5AbsUr1dFHnfbrBTSbgJYugSHacnhA\nCgqn6fspAAKfBFgQwmQaQyDUzc0wdzkQwMMyuUfBFHhaYwEEw+NAjBG0TJjkQV6/6WpAUTJmWpOy\nZtTcykTmfhofttsGESmHB4U7dyfv0L/7HgrGVpslGGXydGk4nPKl3QC0TPjxFMAABIxXk/pwCtcp\nKBDG6xNuBswH3jKAEw24OY24mVQJAnA1jDNbkHaGCBJrG4RgAIiTfrknySrihNtzIs1As8MD0AMs\ndxlHcQxQoDNP/p3b7hprxbtSbDP3GBqqQZKXpJGiLmAFCDNbGAjjFRZQuBLsYYiRkZjZwnhDYB7m\nOOpx4ClWYVpLNYxz4GUEBaLU1sBTRNWYgYN1PU5uEznhQn+6Lh2ZiUcsRuKpp/Z9R8DSK8cAha1y\nRkA5BwBU21bSXUAg0hMX70plEJ/D8j0CwgIOmM6/jH/prg/EBGYCRpoXXDATpn1aMI4DmOLqK8ys\nIA6aiWf1Io477k0bp7acoPl1KJcyilyX1hjFkmdMYme5Wr3W+sW2b0nNAY4ECrfEljYzkoZ4ihaX\nYLFeiSnkjECkmUAQ8xM7AgWAGBZWwJOhUX5f1mQTwDzZFQK6jAOWreAB8GmcgGDEyMNsb4gsIb69\nIzgAaxdmZBAkygO6ypDYJ2IanwoTpm6hb2ENXpDpbackewZ8ebZ4fwnC4bKPIRDGNzLzfyaiFwH4\n7wC+AGGb99dOR8eBiL4TwOsR9uj5V8z81t1GnMluakdH8FTvpC/Wd6oKORNIyhZAIQeIcPp0pj6c\nUhBYAYLc5WVaL8EgjDcE0BAmEgWgGAbGswBOcfLzgNMwLdOemIec7DF9ZYjksM1sjJjMDZZRVDXB\n+j8ULPSRYdQmdN63R2qT/q6Nkh6m8CzCWZHvJKLnA/gtInobgH8O4O3M/H1E9DiAxwF8x24nT2ey\nq82hM3rSNYYdgGBVTmMDWV4x7sAABCh/PKkT81+cN8RJm6vBMs8gMQJhLwYCxum0qRHLoVcTwZi3\nkszfxrHc6j0eVQyEpd7AwiCAtRdC80rkzGLpM2cYg6hjM4oeVcPjam2RHq9LSTznPjwJ4Mnp+hNE\n9D6EQ2Nfg3BIDAA8AeB/AfgOiJOnAfwhEcWTp3+t2M8dvvG7xtAz+Y16dXWhwAryNOVaNSxOR8Ql\nhsaTsDMIWwIPKADC8scjzQfQMAaMNE7/jgE8CPtB2FUWp2FiC8LwGM+OmNnBlKayh+lG44SO5aK0\nGBet/6dmp1iVKdgt1mW3qRitk7zHYNpkUyCiLwDwZQB+HcBj4ui4jyGoF4Dz5Gl56vRDz3mBdwAt\nw13JnpPe1Z4HAJRyGgis6joYggkGymRPvQ95GqdMY0piAQiBLQRgoCG8yXmkZMyTRoFgfACAABpM\nNINDKMeztyIHiDj9EvuCAIkhmzOeuL5S/IOlq6dej7oaspStP1iliX8u8JDiBgUieh6AnwXw7cz8\nFyQeXGZmorZRyFOnn/eil6zr7rBAqhkkzwQEZr29wUDW01SFmC/Ug5X6IPOg1NfuMVMfAlhMuzlh\nyoovdp7sBsCsAozzXtLhYWZeriUwjEyzapGrDZJFANJQuagZgK5KWHEOpfiHkuqRtuFTQ7R+rb5b\n22gVFygQ0T0EQPhJZv65Kfnj8aBZInoxgKem9L6Tp89h6FP72bHtDQAQypFZpggEeVo+gbXruNuS\n8DisWEEOEJnaoD2e4cSpCAhh0tPkoowTMgBEGMg4AsAws4ETwsQNKsTCGoIjgxK1Ivm1EqYQxGIP\nwJpBkMEOcil5H2oBTiU7xbotO77CM8Y9xeN9IISzI9/HzD8ost4C4HUAvm/6fLNIbz95WpFNdoZb\nmPzF+o0gsGqrBgSyvgIQqroggQBrlWFOF2ARVYfq7zSpDwEfeGEPIGAcME6qwnSg9Xx9M3UTp8tN\nMEm4wCGCQDk4OvutJvGARE9YtWXITMoo8NoDGEvd8j+nFTQ8TOHvA/gmAL9DRO+a0v4NAhi8iYhe\nD+AjAF4LYNPJ05uNjQ31e9WAav2dwCBJL6WV2IEGCLJ+rhYgS3MBwcQSYlmOnTMigQAtjCJuFs0s\ntoUTqoHUQqUakVyL7nKqbXklrKAnS82QbVtqRWveXMapeixt2gFZHmlVJTzeh1+F/Xh8tVGn+eTp\nrW/sTe1W2m4BgFBeydAmdi3PYgZ5nnptAIJgATI9T1vZF+a/wkMcjY7TYMJ2kYE1xInMvDCFGPnI\nTBiGMTy8vAQ0BaZAS3QjUg9D8Y0trq0plBsSpasziqVmWH174hZ62cTSbj+r8MhxIhql3NHkr7bV\nCQJqmxUgSNKzT9XIaIHBVE6qDHNebjPIwEPrWxUWn0xAZASMcABN3EyWaQlQmGAiAARjHIOdgSg1\nKEZ1YiChRkx5MhQaKExCcV2cMll1CyS0vrwqhbtM4QfvBQyvHGNnCcr+nJK8ybx99LRl1GVStpFT\nyq4mvYMZbAKEbDylOisQSNI4HYMm1vM5/6CYwWLGDk7VBznxeCo3MiXpcoVknFSz94JDbIMcinWw\nDBvp6nfQanJyNi5Z1z6khor53jL5uM61l+UxmUImTSrRFiZQacPcR1IFDDu/ygxkegkMpu/Jmgcx\nyUtMQaoHc3sJWCgRjLXfLdoYJnYQCMLioiQskY3RMhCBIaoSS18LY+CJIYyI8QpIVIoo0eYwf29k\nDpphUVMxALgZRGzDGpPq7nQwiWX81j3eQpzCbUijPcTNKortFtlDJwgoZZqAQFyXwGAeo/LWTxhE\nsm5hnS9Z2goIYruwZXZLEpBYAad0mr0SE0ggFmWMY7QRLKrE3IDwNqQTeQGHxMaA9SSPt2cZHRPm\noNTNDYf5JKyBRN6GbFvLy/NL5ayywLYFUodQH2ZV1CP5G81qM3sTmm2pdY3dpZV+1T5uARDmMTom\ndN7fih0g/c6UgYsmJZTgdT5HVSJmix9upT4oaoVUJxa1QelalJvVDJFvqRVafzX1AtBVjHwcqzp7\nqhLOci1yKKawkj2YgLOtTaxAKdesPmhg0cIOlGuVORjA5gVbKcQFbEhUiVCYpzRCCFZaysQPMSkn\nD8XcEZAYIOdyUzRkuF7EMkbGdjTJmUPuebDUi9DPvgyiNNatKkdNjgMKDQ/jVhAonjHRCQJquQZW\nkJQpgUH8rrztc5VAA4SE2ViMQKx1qNsRwl+wHSyKwaJCYPZEzN+BGRhofvvHxuSPInURzHYGQFEn\ngGQLNxnfEL9H0dSKXDwAAaW+tflLlBaQkO1J6d0sxiuHUB+8UlUHgOqDXGQEG1jBrQKCInkfqs1A\nq6/Vy9su5BVFPqtGZRagkE+YqEok6sX0qXkUcnVCUmuNyudTKVctLK+Fll9KL3kKSmpGqR+tzF5q\nxHGYQiZnVQkq9ZvVA6VM0c6QT3ytbvLdYAhTWlFlyMqU6ln9m/dhSVQLJtIAJizbtQE0+x4WNSIU\nFUwj6WjNGOJbXlMngIU1yCXXJdYQb620N0POHGR+3pbKKgwVAygziLw9rb/WciU5DCjsAQKhnUKh\nFhAwyhcne55fYwUiXU/L2EE+kZFN4BwQamPJwCPUEQ+qAQ7m75V4HiDAAHN4M6aAphkYaHpjK8AQ\nbz9MmAlG4uTMJ202lBwcNHtDHgotlZeSLWEPgAhj6geJvN287VK5mtwf6kPpzTWJ6TGo1Dff+Mpk\nPxcgaG1qgLAuk+ark1hjCVDSauOqieNlxFqZwgMrVQepSkjKrakTsXxIT79rZfL6JbXC6muLejHn\nbVAzam23yGGYQiJbGUGhjU2sQCnnURNW5bT0Chh4GUKSlzMHrK/TT1YBpPr/yBmCvI7qgFAhAh8Q\nbMFUI5YGZ4YwMYhkgijqRCwf2EAIk1bVB1HeeqN7mAOAXdhDGJ/IyyCqh0FofZTkOKDgYAK99VuA\nQC1fA4xCeS8YhLQ6IGgTesUAciDIxWIdjnJuERQ/fA+NcFQLohpBIvJxBQyxYxsYmtUJAQ6J+iBv\ntaBWyD5i/p4AscorqBh5P7KvXFoYxOHVh6JaAFRVA696YJZXmMFdAcKqLY0BKOPR1AatbbWvBiHt\nx66pDM4XmKZKxGv5ua6X5ksPhak+ONWKPF/7bo2ttlbCktqaB4+aUZPjMIVJtjCCUL+tXpeKkJcp\nlW8Fg+m6OmFLgKABhAVe8yfr94KK5IwgViL5OZWJ6oRQI2btAmu2EJqSLAHi2skYhBoh1YmQZ7AG\nYUi0WIP8qSyPhUzTmINVtpQ+5xdUjLy/vM+aHIYpVBkBcAxAKDCTVXmLNcxpSqLBYJI8DRCM+sXx\n54BTGYdLrGcveUPLa19HeSh0/LQYQ26A1BjFFtYg+7DyrbTWlZbVsOgKe4h9euUYTKFxMkcp3ucG\nIFDLVb67wGCehLQum9cpMAR1DHmdDCzy+yn9DqYaUXrb5KEFWrq4nl2UmN70itExDoEhjYsKOxBa\nv2QEsbs8TmG+nQ2sYc6T7cGITzDe+lvYQ54H1O0PXjkGKERxgNluQKCUdwFG6e1fAYOQngGCVqcA\nCGoZBVDU/vM8cV20J3SzBlLUBnk9TWg5fmF0JNFODRikejEPW0xaS52I5aKHIoyBEuNgrgLkhsM5\nP/u5zgkQtTygrmJYchj1wcMKipPbeNubE70GCFqbTkCw61C9/waGkE/8BoaYtdf3RpmFjc+G+nIV\npWyWBWqlakf8pOQzqV+g85o6MnIa15Bv4mK1W1MrtDJWmjVubRzePGDNIkpyLKYgpPqAb2QFZvnK\nm3IXdiCudfUiK6uwgqRMpc1VPc/zUWIOHsmNjHPnPKfHaIXcdRnVCIg3NM3XYVCzKzNjDJoqEb7y\nfDtSnSCljz1ZQ+xTKyPT8nSLOZTq1PK8cihQ6AWCYl2DQVTLeZiElucEA62eRfmLKkOWb41FG0No\nR3gdCkDgZSHEk00gqZy1m9sfpskc4xZkujw4hkIHmZpQVyUsOwOAVbBTvK7ZGkKebyWlR60opd82\nQFTVByJ6CRH9MhG9l4h+l4i+bUr/biL6KBG9a/p7tajznUT0QSJ6PxF9rWcgVVtB4U3vVRHMfjYA\nQrO09L9hciblK+xAtplct3XVJgz1ZqIaoYZEW01VVAnLM6GVya8tD0XI0/sx85V+Sx4HTWorKnvy\nctly6jQA/BAzf78svNup0z2swKjnZhGNYNDEELSJV0pTWIFWpsQS8vF7AOKsMr+9F+Oi9ELkDCKJ\ndASwjl8I7dVUCZVJIGMMGVPwsAZAURs61AqtnJUm+5D95HW0el6pMgVmfpKZ3zldfwJAPHXakvnU\naWb+QwDx1Om6VB7aKitQJoE6sfOy2fdVvWxiJvlJnoi1kG0WJv+qrelzTi8AQjK+/B6zsesAwUl6\nKwuxZBXZOL39LWHJGjSjo/YWV9LMmAUlnTnsFZl7J5I+suucNeQxDHkbpfzpVtVYBy1K0tozYb4X\ng230LJJq8j5kp04DwLcS0buJ6MeJ6IVT2ucB+CNRzTx1mojeQUTvePbpT1WZQYuKEOuo5StpTaqC\n57euMIRVOU9/63lhswKP7MEe8qe79pKKM6L0wEo1guUW8XLyrdOSJioUXw5VTqASMOypUsj+83Kt\nYdC9qkUublDIT50G8CMAvhDAywE8CeAH3L0CYOY3MvMrmPkVV488d52fv5GTwcBUE4rsIP/uZAdm\n/pxXZggqG4BC5w1WoPVbUhtMMViDVq6LeNaAwWqUKWULKD/gsl4ou9SJJ1znkzo9+VqZ+KJZCQw5\nu1jK6O7LWN/DGvJ1FV7mYLUh+2qJYMyl+9RpZv64yP9RAL8wfe07dRqVh7XCJFzlDSAplSnl17wL\nSf0aYygBQgE0TMDIy6zuQznXIR9TPt5VhlJeS8vzQTBjI5hCIQJkpGOoldoX5uGqNoXlZnKvxNr2\ngCTQieCzNQBr92Wat7Sh5Wtl5p8I6c/o8TK02B5K4vE+EJRTp6fj56N8PYD3TNdvAfANRPQwEb0U\nzlOnewChyCQcabcGCKU2RZ0iQ8jrF34TdVy9srWNAjtQX41nFMvGkF+3Sr43ZKtKoZUB7J+mJ4Bp\n77UP1qnT30hEL0cY+4cBfEvovP/U6UQ8D72nzm2BwarcOr2ZIciyLW9wq8zMGrLHzWQT9fbdwiV2\nEDqb91qIadl19EZAeA/m+AUgCYVe2EJoOw9SkowBwMozAbQzhlAOcz9p/jbWEMeTl83LW+20yJZT\np3+xUKf51GnA6CW2aeU50z2Mog4WjexAS98ACHLiVkGH1t/PIWE5Q2HCS4mqQRyTVkRMfExtU1wz\nQaIdAQx5/dhBfqw9oKsSMV/mAWmXeTQjGZO7Fxy0MrJcHA+Ue7ZA4GzBS7cmewGCMgHOCgiWWECR\n5ZvqQIEhePsticoM9hJGcEt6nsMd1AfNTSnTtUAoT14+vJLxcK0ylMfs9QbsoVa09AccJcy5hRo3\npHeBQVZGPVVaLVdJL73ZS2NR2nDbKrSxKKpD8fu5JDc4iu8ltiCNjgCQhEETLyxCMTzK74CeB6zV\nA0udkG3k10DKGkqMwLsRy1a1wivHYQqTFI2HDcxAffN72ME88WjNDhoAIRmDg+qr4/BMfq2f/P6z\n77WXhjDc7yO1Z1PmGwMtGco4L2MYFFPXJal5eV8l1uBxXYay2XjZEdiktqMbEi177a0EL51bbpMd\n1NQFMy+brBZDKLVTfNNbwOX9v7aMJ++j2G4nz9eqbVEZ5okdv1OxyVw98KoSPd6ILepEKKNMeCXN\n9DA0tm3JIUCBUXhrO9P3YAchX2EIyvWKCVgMIZvUFqtI2jfSi0DW+gw7fqtqnR4p6UzJpBT5uQsx\nBwbZXMG+4HFJasCQ1E2HWXQ3lgKetPzYXgtrsEKnt+DuIUAhEevhLNBpLxisJmQNDLTJDZsdmKxB\nm+QWIORtVQCBLeARaUtbvBpjxwvRlpanUZYVE1++IZJJX2o3YwuWGiHz9lAlYrFaBKPGGvYCB6ts\nPsYWORYoWA+oke5RFdRyG96ITeqCKFNUGZx9e8SnCmzr4y7FZAvWG1sDGEf7HmAQwwBQpuheW0Mu\nXpWi1H8rMBwHFJzsIL7xu9iB0o+5bkH5XmIIJZWh1H+SZrCHIksoSTLejneGF8iMDPUcCK+w8sMb\nakRSTUs21IbaxNfWS8TrVgOk1k8ov401tBgjvXIcUJCiTHCg8Gb1gkHynBVsB8jLGpM+H1MBNGTa\nqr0sT/2OcjnVC9EhO4QN6A2ZHgaljGVbsLoqsIX1pFzy4vda2LOmaqzKZLdSUydaPRSxTYs5eG0O\nHjkWKBTAoFtVUMoVYw+MdjXpAYSefqpSaFt74baMYxebg+kaqHSkpk0flhphdemk4TVgKJXJpbbm\nweOhaPI+NKoVlhwHFIxxe9mBWlZlEQVAUFhASWUw28jLlNJEn90sQSmTluN1fuk52QoEsjuTFVTY\ngkzKmYDyZk/atSZfQW1wGRcbDZCADxhKHgqtTGzXyxq0cZTkOKAgJP5ffZPcmLwKGHiDkWK7WtmS\nKtEMCFp/FiBYbVTKr+pZdYogsZtS4Zc8DqHofYBaxmNfCNdLmleV0Mp47Axa2bytUOc84OCRw4FC\nC8X1qxQFdqB8d08Wo40iIJTG4Wz/XGXZSL812aKrrGwY5UkX0ipNNgJDbUit6kSoUx6X1bZnfJYc\nBhSKdoNOQFDPpywAQkldWPVZYwi1ehbz0NrI061xqb/JHbzlW6QU8ajmlVUINaAJWVllgjerCYbH\noMQYAB0YPEbI22QNhwEFVYwJ5gWEantFyrzuV82zAKEVoHtekoXfYctL906kx8tgeDg8Hguv4bEk\n7tiH7HuPR8Aeg581eOWYoNDKDmqAkJfJvpcYQpE9OABBtT8UJq6bJVTYQ0jjepnbEsdb365bHnhx\nMsc3uFa2wb5gMoEGG0OLncHLGOyyfe5I4GigUAADdbLU1AWtvRoDKKkCWwDB6F+mWf/D6v/W8793\nANLuUlINvHWtt76lRhhGx6Spin1hL1VC66vHztACDrnc36skWyaEChyOGy8BQq3fXkCotd+gwlQn\nf++EPzNAbIpszMRF2StGx9DO+sdqOZUqaaNSdyswWGle1qD1UZJjgMIGdhDKVtSFrA+zbZFvjq8V\nEBT1oUVtUKWlHUXOqk44XIfVfNNmkNF9iy3EbF5PVnUIyuROVQY9Pe+z5JXYaoC00lpYg1eOAQpC\nWl+X52EAAA0uSURBVMHA5V2oTfhzAwLK+S7GYtF/q25uT6jI/FztBRRuNSEDgF6VwzA6Jv047Ash\nXSQbP2AJGCwAWZXNhmmtcSj1vdT1eSg8cihQaBn/Zu+Cku8FBJfU6raoDa393UY9r+zlETVDpeuT\nqNhsQ/n1hi06MLT20atOeMfe+pt4zn14hIh+g4h+ezp1+j9M6S8iorcR0QemzxeKOvucOq0wBJMd\nVABhxUA6ASFpp8YylLpmmiIWSyj10VX/tsWh87tVCK35ktExYwtpeVLb77Ed1ACkx87gYQ17qBMe\npvBpAP+Imb8U4Yi4VxHRVwB4HMDbmfllAN4+fc9PnX4VgB8mopN7RMB220HDhFfBokNlUIGiMvlL\nIFJSD6y8Jb3z9WypKLVq54qPsgZQsS2kZZGWXWdN+dqE0/NrtgNZzls2jsejTmyxNXjEc+o0M/Mn\np6/3pj9GOF36iSn9CQBfN133nzoN+G0H7rr1MmbeOd6sG9WGbvuRA2w2ywbjVqiffVr5nqZKk7w2\nTgcAtEz2Ut4Wg2CtbaAPGFw2BSI6TadDPQXgbcz86wAeY+YnpyIfA/DYdN1+6vRffcrPDoB+daHE\nILxsoqQyeBhGIa2YZ41vlX6AgCXlQZzdkrWJ72xvzmpgC2laQY2Q7ZWG1cAYrDz1e9ZPy9Zse7AG\nFygw8w0zvxzhsNhXEtHfyfJz5uNpczl1+lHt1OkG20GDuhDzWwDBkiZAEON0TXqZ5qDzTS+cuwCK\nmpQmtZZe/DGmj4pKoDVtRTzaLkmjbvZdUyVaIiCBdnDolSbvAzP/PwC/jGAr+Hg8ZHb6fGoq1n3q\n9NyPR1UAfOpCaxtOet/LEGriesPXQNArRwSHXPKbK8ZBFCZHgS3Umm2RrepBjTUA/n0be4HB4334\nbCJ6wXT9KICvAfB7CKdLv24q9joAb56uu06dBgq2A4e6EOrX63WpDEpeTZrqtjIIK886/aky3nPZ\nCXNpimz0BEEpnoJi/YLRMSlmqBH2yshtqkSrOgG0AUMrOHiOjXsxgCcmD8IA4E3M/AtE9GsA3kRE\nrwfwEQCvDYPoO3XaZAc9YKCUqxocSxPXwxCMsRXHYNUp3aNjoleM9i7Z7MkAwtNM4rO3vvVdZjHC\nYbFMqzEzp8e5rcoDyA+sTeqINuWhtHk57cBaqxxgt6N+nz7l7avHyU1p6/v1/wM8p06/G8CXKel/\nBuCrjTp9p05L2UJtNwBC9zi8wNEhOxqojyXWJF+BwXqiu/Ja+lSyvcBgtlWZ6K31HcPv7kvKoSIa\nZ7Henspk9xgVi203qBfNdgSn2uDJV7/LehWvg8ke9OT6WLaIMzqxqw3ZToHSh7TOvrO6vbEJNeOj\nOe7s+95bsh0LFJSJP6dnsrfKsCqfTdjdAUHrV97/HuxFK2vVPTob0QyGK/1cqyfzNbAoVim6KcuT\nv5Rnl7XK1wKdgP22ZDsWKGjiAQSDRRTbcb5Ji+PZUcWxxDWunQ56KabflrR4G6w6VlqtXcMb4Y1f\n2BK0VAMGr5RYg1eOAwrOiV6d7NgBEDrsA2dVG2Raoe5d2h5MXOqNTvR4H5K0cPNuT4TMNtQIDzCU\n10Lk32+HMQDbtmQ7DihI2aAubFYZvLR/T7VBGVeS3/oi9KggdyQrt2RTZKP2QyhJeZRjkif/QUqd\nUvcOYCgFN63aUL+Xy5tpyngfjBOigG5AcLXV8Gb1MJIuuYOJerZ7ObfsGUDRCT7Fap2MYZ1/PmDo\nkeOAgqIuAH77wZ4qgzo2pf+SAbIJoJU2/GqFEbBU6+fcsgpLbqnrYBMVm4CaJkKVtXaa1Yi8qwbD\nYw9j8HomvN4JS44BCsb/stt+sBEQPJ6GlXgAQStTa9/JnHpUyNW8uisG4TImxk9lQpdUCKu5vYCh\nYfJvVSWsOnuzhmOAQqf4LPOVOg0qhVmucRxNtoZeuV9UhNuQc1tgG9rfqkpoZay0XjkcKJgMwTN5\n92YIrRPfGpdHavU0taJlH8bbUitafoAWhlDrQ2ES7riFVbreV84WvB6JmqtyD8ag1oPvJ87lUKCw\nycOwJyAUxtBlR7DKOIHAZUswynWpFS11mmwFZxjEBttCj/QaHkuRjHm+/t0HNpbbskUOAwrdHoYe\n8U7Gre10SMkN+aBI9zkQHZN5b9tCsX51LG337dkPwttmy093GFBYyRajotdwiPKbP//e5MKs1Veu\ne2wNJTC9VaDdKq0TvqRCeOrspEaU+tDe7knRXg9DZz2vHBMUvCpDT70sv9hmg0rRVKZHhTiXbLVH\neKSXrscfrSf6sWZb2DKerMuSfSHvv0ct8ACDVe/BXRC1Z/udZc9twC71vc7f/pTvGRd0Z7LVYNLK\nFvKm2ntf97WhjvvU68b+jgMKxrh3URmU9ov5NfXjDCxhs5uyaJSs1L0DWW3kWpNCTIKZfw4k99oX\nKh6E3YKVHPWsupYcBxQU6bZHbQGELf3sLeoEX56C8/d/n/KJjmCmYnO9bKERGHxjOT/CHxIUTFvA\nuRlCrbz3/9H4f6v14WYJLXUcbXTLng9uj0HRqK8ulCqpEHt4IxwxBzKvZ92DVa8XQA4HCiYYON7u\nZwWEWl+eNpTrXgNk1/+7kbFvFveahUJetQ99gvfWc0+kzOhYXB9RAZHdFkQ5g5xqcihQ2PKC2ZtV\n7aFidHscLKnR+QPaDnaVcwVKOdtZTbrOh65HldliJ7hvDY1NDMEDxBtYQqlszc25t2zu634DilbD\no6usokJk6V1soTCkVjViSwRjrZ5V15Itp05/NxF9lIjeNf29WtRpOnVa/f86J7rHO7FXPEJVLGAx\nrndhEmeyc9yGuCMbW1wz59CNMrZwLvtCbwTjFnelJp5zH+Kp058konsAfpWI/seU90PM/P2ycHbq\n9OcC+CUi+iLP2Q9LI+6SVbkVO8Ity2b3ZZSD3E+3MMr3IPOZdPVLpovyzNkW6bW+jGLrdoxxWOWd\nZfQ0uLaiz2XLqdOW7H7q9G4qQ2v5vVlCo6TtdED/VjvIHrLFiFhqb5VusAWPe3JHo6PWZW1jltYY\nBqvNvRjDllOnAeBbiejdRPTjRPTCKa351Ombv/rUlKh1vk7qVhlaWEALwHgBoVd1sGITlOvdJvVt\nxyh0rlto68OwLbjri8uKGlGNX+iwMayaORMwbDl1+kcAfCGAlwN4EsAPtHQsT50+Pfpc91ttv4fe\nX3SXPrfYG3aU+zQc6c5kr2ChPQKVbgsYuk+dZuaPT2AxAvhRLCrC5lOnAWxSGTbbEfZQG25T7jfb\nQM8CJ09ZjwrR6Imoja2FLWwJbCqV2Xsnpu5Tp+Mx9JN8PYD3TNfdp04vna6TPIDgaqsRELxqQzFv\nC3B4VIdqGxvzvd3sOdnvWkq2hSJo7AsMW4KUeoFhy6nT/42IXo5w3x8G8C1hIH2nTpfECwibA45u\n+417y6rDEYWYwF77RcVyv1kK7Ze8Aisr/8Zxlk62Nvs0y7UfNLvl1OlvKtTpO3V6A0NwAUIDC2jx\nXOwR/biHbPV63FdglLgb5bWYkI4yxQldckNmeSX331Y3pVrH6HMPYDhMROOugOBofy87QjXPY2Bs\nlT3auWsAKK2B2NKGt3yHJ+K21Yi7UiWId9uWpl+I6E8AfArAn971WG5RPguX+32Q5Yj3+zeZ+bNr\nhQ4BCgBARO9g5lfc9ThuSy73+2DL/Xy/x1EfLnKRixxCLqBwkYtcJJEjgcIb73oAtyyX+32w5b69\n38PYFC5ykYscQ47EFC5ykYscQO4cFIjoVdNmLB8kosfvejx7yLRq9Ckieo9IexERvY2IPjB9vlDk\nNW1KczQhopcQ0S8T0XunjXi+bUp/IO+5sPHQg3G/zHxnfwBOAP4AYbXlQwB+G8AX3+WYdrqvfwDg\nywG8R6T9JwCPT9ePA/iP0/UXT/f9MICXTr/H6a7vofF+Xwzgy6fr5wP4/em+Hsh7Rgj9et50fQ/A\nrwP4igflfu+aKbwSwAeZ+UPM/AyAn0HYpOW+Fmb+FQB/niW/BsAT0/UTAL5OpPdvSnMAYeYnmfmd\n0/UnALwPYQ+NB/KeOYi28dADcb93DQquDVkeEHmMmZ+crj8G4LHp+oH6DYjoCxDWyvw6HuB7NjYe\neiDu965B4a+lcOCUD5zbh4ieB+BnAXw7M/+FzHvQ7pn1jYdk/n17v3cNCvtsyHJ/yMfjHhTT51NT\n+gPxG0yb+v4sgJ9k5p+bkh/oewbSjYfwgNzvXYPCbwJ4GRG9lIgeQtgF+i13PKZzyVsAvG66fh2A\nN4v0bZvS3LEQEQH4MQDvY+YfFFkP5D1bGw/hQbnfu7Z0Ang1grX6DwB8112PZ6d7+mmEfSuvEfTH\n1wP4GwDeDuADAH4JwItE+e+a7v/9AP7JXY+/436/EoEqvxvAu6a/Vz+o9wzgSwD8n+l+3wPg303p\nD8T9XiIaL3KRiyRy1+rDRS5ykYPJBRQucpGLJHIBhYtc5CKJXEDhIhe5SCIXULjIRS6SyAUULnKR\niyRyAYWLXOQiiVxA4SIXuUgi/x/eJGSpNSyrFQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10bc1d68>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lat, lon = get_latlon()\n",
    "\n",
    "plt.imshow(lon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.75 0.75 90.0 0.0\n"
     ]
    }
   ],
   "source": [
    "rlat = precip.coords['lat'].data\n",
    "rlon = precip.coords['lon'].data\n",
    "drlat = rlat[1] - rlat[0]\n",
    "drlon = rlon[1] - rlon[0]\n",
    "rlat0 = rlat[0]\n",
    "rlon0 = rlon[0]\n",
    "print drlat, drlon, rlat0, rlon0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13.333333333333334"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "j = np.floor((lat.flatten - rlat0)/drlat)\n",
    "i = np.floor((lon.flatten - rlon0)/drlon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def regridll(ilat, ilon, olat, olon, test=False):\n",
    "    '''\n",
    "    ilat - scalar or array of input latitudes\n",
    "    ilon - scalar or array of input latitudes\n",
    "    olat - array of output latitude\n",
    "    olon - array of output longitude\n",
    "    '''\n",
    "    \n",
    "    dolat = olat[1] - olat[0]\n",
    "    dolon = olon[1] - olon[0]\n",
    "    olat0 = olat[0]\n",
    "    olon0 = olon[0]\n",
    "\n",
    "    j = np.floor((ilat.flatten() - olat0)/dolat).astype(int)\n",
    "    i = np.floor((ilon.flatten() - olon0)/dolon).astype(int)\n",
    "    \n",
    "    if test:\n",
    "        print 'Input -> Output'\n",
    "        print 'ilat={:9.4f}, j={}, olat={:9.4f}'.format(ilat, j[0], olat[j])\n",
    "        print 'ilon={:9.4f}, i={}, olon={:9.4f}'.format(ilon, i[0], olon[i])\n",
    "        \n",
    "    return i, j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "i, j = regridll(lat[0,0], lon[0,0], precip.coords['lat'].data, precip.coords['lon'].data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-135.0, 29.896942, -180, 80, 135.0, 30.0)"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lon[0,0], lat[0,0], i[0], j[0], precip.coords['lon'].data[180], precip.coords['lat'].data[80]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "300"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ilon = 225.\n",
    "(np.abs(precip.coords['lon'].data -ilon)).argmin()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "135.0"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "precip.coords['lon'].data[180]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "alon = np.array([0., 45., 90., 135., 180., 225., 270., 315.])\n",
    "alat = np.array([-90., -45., 0., 45., 90.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getij(ilon, ilat, olon, olat):\n",
    "    \n",
    "    # Make sure everything is a numpy array\n",
    "    iilon = np.array(ilon)\n",
    "    iilat = np.array(ilat)\n",
    "    oolon = np.array(olon)\n",
    "    oolat = np.array(olat)\n",
    "    \n",
    "    # Convert longitude arrays to span 0. to 360.\n",
    "    iilon = np.where(iilon < 0., 360.+iilon, iilon)\n",
    "    oolon = np.where(oolon < 0., 360.+oolon, oolon)\n",
    "    \n",
    "    i = np.abs(oolon - iilon).argmin()\n",
    "    j = np.abs(oolat - iilat).argmin()\n",
    "    \n",
    "    return i, j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 45.0\n"
     ]
    }
   ],
   "source": [
    "i, j = getij(1., 27., alon, alat)\n",
    "print alon[i], alat[j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 355.  310.  265.  220.  175.  130.   85.   40.]\n"
     ]
    }
   ],
   "source": [
    "print np.abs(alon - 355.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
